Changes in annual transcriptome dynamics of a clone of Japanese cedar (Cryptomeria japonica D. Don) planted under different climate conditions

Environmental responses are critical for plant growth and survival under different climate conditions. To elucidate the underlying biological mechanisms of environmental responses in Japanese cedar (Cryptomeria japonica D. Don), the annual transcriptome dynamics of common clonal trees (Godai1) planted at three different climate sites (Yamagata, Ibaraki, and Kumamoto Prefectures) were analyzed using microarrays. Both principal component analysis (PCA) and hierarchical clustering of the microarray data indicated the transition to dormant transcriptome status occurred earlier and the transition to active growth status later in the colder region. Interestingly, PCA also indicated that the transcriptomes of trees grown under three different conditions were similar during the growth period (June to September), whereas the transcriptomes differed between sites during the dormant period (January to March). In between-site comparisons, analyses of the annual expression profiles of genes for sites ‘Yamagata vs. Kumamoto’, ‘Yamagata vs. Ibaraki’, and ‘Ibaraki vs. Kumamoto’ identified 1,473, 1,137, and 925 targets exhibiting significantly different expression patterns, respectively. The total of 2,505 targets that exhibited significantly different expression patterns in all three comparisons may play important roles in enabling cuttings to adapt to local environmental conditions. Partial least-squares regression analysis and Pearson correlation coefficient analysis revealed that air temperature and day length were the dominant factors controlling the expression levels of these targets. GO and Pfam enrichment analyses indicated that these targets include genes that may contribute to environmental adaptation, such as genes related to stress and abiotic stimulus responses. This study provided fundamental information regarding transcripts that may play an important role in adaptation to environmental conditions at different planting sites.


Introduction
Based on analyses of tree-height and tree-ring data, many studies have reported that tree growth is affected by environmental conditions [1][2][3][4]. As Japanese cedar (Cryptomeria japonica D. Don, also known as 'Sugi' in Japan) is a major forestry species in Japan, accounting for 44% of artificial forests [5], and is widely planted throughout the Japanese archipelago, understanding the effects of environmental conditions on this species is very important. The relationship between environmental conditions and phenotypic traits in Japanese cedar has been examined in several studies. A positive correlation between annual ring growth and temperature during February and March was reported by Takahashi [6], and significant correlations between tree age-height and several climatic variables such as warmth, solar radiation, precipitation, and snow depth were reported by Nishizono et al. [7]. Quantitative trait locus (QTL) analysis of several phenotypic traits in three replicated common garden experiments under three different climate conditions identified an average of 53 QTLs [8]. However, only two of these QTLs affected the same traits across all three environments [8]. Moreover, one QTL associated with growth response to drought stress was identified by QTL analysis of trees growing in two contrasting environments [9]. The small number of QTLs common between regions and a QTL with a large contribution to overall climate sensitivity would indicate that Japanese cedar is highly sensitive to environmental differences. Although these previous reports suggested that environmental conditions affect the phenotypic traits of Japanese cedar, there are no published reports of studies examining the underlying biological mechanisms.
Japanese cedar grows continuously until environmental factors (temperature, day length, etc.) or internal factors (aging, etc.) cause growth to cease (indeterminate growth), whereas the degree of annual growth of many other coniferous species, such as those of the genera Picea and Pinus, is regulated endogenously (determinate growth) [10][11][12]. As such, environmental conditions may exert a greater influence on the timing of growth cessation in Japanese cedar than in determinate growth species. Results from studies using growth chambers indicated that a short photoperiod and low temperature suppressed height growth of Japanese cedar [13,14]. Prior to growth cessation in autumn, Japanese cedar trees accumulate starch and break it down into sugar in all parts of seedlings (upper, middle, and lower layer shoots, and roots) to enhance cold and frost tolerance as a means of ensuring survival under harsh winter conditions, and this process also contributes energy for bud breaking and shoot growth in the following spring [15,16]. Soluble nitrogen and total free amino acids in the lateral shoots, stem, and roots increase at the time of growth initiation, but the levels decrease beginning in March in the roots and in April in lateral shoots expanded in the previous year, and in stems [17]. Since new lateral shoots exhibit high accumulation of soluble nitrogen and total free amino acids in May, stored nitrogenous compounds may be transported to new lateral shoots for growth [17].
A few studies of evergreen coniferous species have used time-series transcriptome analyses to investigate transcripts of needles and shoots, the key perennial organs that sense changes in environmental conditions, and these analyses revealed dramatic annual dynamics [14,[18][19][20]. The annual transcriptome dynamics of Japanese cedar were also clearly demonstrated in our previous report [14]. Principal component analysis of microarray data demonstrated the seasonal cycle of the transcriptome and explained the seasonal phenomena of Japanese cedar. For example, the expression of growth-related genes was up-regulated during the active growth period, and the expression of genes associated with starch metabolic processes that may contribute to cold tolerance was up-regulated during the dormant period. As the respective effects of day length and temperature interact to control this annual transcriptome dynamics [14], the dynamic pattern may change under different environmental conditions, leading in turn to physiological and phenological changes that help trees adapt and survive in different planting sites.
To investigate the changes in annual transcriptome dynamics of Japanese cedar under different environmental conditions, we planted specimens of a single clone at three sites covering a wide range of environmental conditions of natural distribution and analyzed the transcriptome. In this study, microarrays were used for transcriptome analysis. Although microarrays cannot detect novel sequences and splice variants as RNA-seq can, it is a high-throughput, cost-effective, and highly sensitive technology [21] suitable for the purpose of this study. To the best of our knowledge, this is the first report of environmental responses investigated by analyzing annual transcriptome dynamics in multiple sites in a tree species. Our data identified differentially expressed genes that may be related to molecular mechanisms controlling physiological and/or phenological differences in a clone of Japanese cedar between planting sites.

Test sites and plant materials
One-year-old rooted cuttings of a Japanese cedar plus-tree clone, Godai1, were planted from March to April 2013 at nurseries in three locations, Yamagata (Higashine City, Yamagata Prefecture, Japan [38˚23'53"N, 140˚22'47"E]), Ibaraki (Hitachi City, Ibaraki Prefecture, Japan [364 1'28"N, 140˚41'21"E]), and Kumamoto (Koshi City, Kumamoto Prefecture, Japan [325 2'53"N, 130˚44'5"E]) ( Fig 1A). The clone 'Godai1' was selected from Kimitsu, Chiba (351 3'12.0"N 140˚07'48.0"E) as a first-generation plus tree. We have been studying 'Godai1' as a model clone of Japanese cedar [14], as it shows average growth and high rooting ability for propagation [22]. A total of 10, 6, and 8 rooted cuttings were planted at more than 25-cm intervals at the respective sites. No obstacles hindered sunlight from reaching the rooted cuttings at the nurseries at the three sites. Data regarding air temperature and day length at the three sites are shown in Fig 1B and 1C. Yamagata is located in northern Japan, which is a cold region with snow cover in winter; Ibaraki is located along the coast of the Pacific Ocean and has a mild climate; and Kumamoto is located in southern Japan, which has a warm climate.
Throughout the year of the study, Kumamoto exhibited the highest mean air temperature and Yamagata the lowest mean air temperature among the three sites. The mean air temperature in Ibaraki was intermediate between that of Yamagata and Kumamoto. In Ibaraki, the temperature in summer (June to September) was similar to that in Yamagata, and the temperature in winter (January to March) was similar to that in Kumamoto. A 10-cm-long apex portion of sunny upper branches (S1 Fig) was randomly collected each month from two of the planted trees at 10:00-11:00 AM from February 2014 to February 2015 (S1 Table). All sampling at the three sites was performed within 7 days. Meteorological data for the sampling days are shown in S2 Table. The height of all cuttings at the three sites was also measured at each sampling time point and other times, and growth rate was calculated by dividing the height growth differential from April 2014 by the initial height in April 2014. The significance of differences in annual growth rate (growth rate from April to December 2014) among the three sites was tested using one-way analysis of variance (ANOVA).

RNA extraction and microarray gene expression analysis
All collected samples (3 sites × 13 time points × 2 ramets) were immediately frozen in liquid nitrogen and stored at −80˚C until analysis. Total RNA was extracted from samples as reported by Gehrig et al. [23] using an RNeasy Plant Mini kit (Qiagen, Hilden, Germany), and DNase digestion was performed on-column using an RNase-free DNase set (Qiagen). A NanoDrop 1000 spectrophotometer (Thermo Scientific, Waltham, MA, USA) was used to accurately measure RNA concentration. RNA integrity was assessed using an Agilent 2100 bioanalyzer (Agilent Technologies, Mississauga, ON, Canada).
The microarray probes were designed using e-array (Agilent Technologies) with the default Base Composition Methodology based on isotig sequences determined from next-generation sequencing data collected from analyses of various organs (cambium region, sapwood and heartwood, tree-top, shoot, male strobili, roots of seedling) in several developmental stages Sampling sites (A), mean air temperature (B), day length (C), and annual height growth of Japanese cedar (D) at the three study sites. We edited the map of Japan obtained from the GSI Maps Vector provided by the Geospatial Information Authority of Japan (https://maps.gsi.go.jp/vector/). https://doi.org/10.1371/journal.pone.0277797.g001 and seasons using a Roche 454 Genome Sequencer [24]. A SurePrint G3 Gene Expression Custom 8×60K Array (Agilent Technologies) consisting of three probe sets corresponding to 19,304 sequences was used for microarray analyses (GenBank accession number, GPL21366) [25]. Gene annotations shown represent the top-scoring BLASTX hits for each sequence's predicted protein product as a query against the TAIR Arabidopsis protein database TAIR10-pep-20101214, with a threshold e-value of 10 −5 , using the CLC Genomic Workbench software package, version 4.1.1 (CLC bio, Aarhus, Denmark).
A total of 200 μg of total RNA was transcribed to double-stranded (ds) cDNA, and the cRNA was amplified and labeled using a Low-Input Quick-Amp Labeling kit (Agilent Technologies). The cRNA was hybridized to the custom gene expression microarray for 17 h at 65˚C and washed using a Gene Expression Hybridization kit (Agilent Technologies) according to the manufacturer's instructions. The resulting slides were scanned using a SureScan Microarray Scanner G4900DA (Agilent Technologies), and scan data were compiled using Agilent Feature Extraction Software 11.5.1.1 (Agilent Technologies).

Microarray data analysis
Raw microarray data were normalized using a 75th-percentile shift and analyzed using Gene-Spring software, version 13.1 (Agilent Technologies). To enable direct comparisons of transcript profiles, median log 2 -transformed ratios for each time point were normalized to baseline (normalized intensity value). A total of 13,385 targets with a raw signal value exceeding 100 in each pair of samples under at least one of the 39 conditions (3 sites × 13 time points) were selected for further analysis. To assess the overall trend in annual transcriptome dynamics at the three sites, principal component analysis (PCA) and hierarchical clustering were performed using the 13,385 targets for the 78 samples. PCA was carried out using GeneSpring software. Hierarchical clustering was carried out using the R package "pvclust" (distance: correlation, cluster methods: average, replications: 100) [26].
To detect differences between the sites and identify statistically significant differences in profiles in the time course transcriptome data, the microarray data were analyzed using the R package 'maSigPro' [27]. Three between-site comparisons ('Yamagata vs. Kumamoto', 'Yamagata vs. Ibaraki', and 'Ibaraki vs. Kumamoto') were performed using 13 time points of normalized intensity values (cubic regression, false-discovery rate �10 −3 , R 2 �0.7). Although a limited number of replicates from each time point were examined (2 replicates), we performed the regression analysis using a sufficient number of time series samples to detect significantly differentially expressed targets.
To identify environmental factors that affect the expression of all 2,505 of the differentially expressed targets among the three comparisons, partial least squares regression (PLSR) analysis was used. In the PLSR analysis, target expression level was used as the response variable, and the meteorological data were used as the predictor variables (S2 Table). As these predictor variables are mutually correlated, we used PLSR to account for the problem of multicollinearity. Meteorological data for the observation sites (Yamagata, observation sites in Higashinecity and Yamagata-city; Ibaraki, observation sites in Hitachi-city; and Kumamoto, observation sites in Kumamoto-city) required to calculate air temperature, precipitation, and sunlight parameters were downloaded from the Japan Meteorological Agency website (https://www. jma.go.jp/jma/indexe.html). To determine the effects of short-term and long-term temperature on the targets, seven parameters of air temperature were calculated across the sampling sites and dates, as follows: air temperature at sampling time ('temp-0h'); minimum air temperature on the sampling day ('temp-min'); maximum air temperature on the day before the sampling day ('temp-max'); and mean air temperature of the previous 7, 30, 60, and 90 days ('temp-7d', 'temp-30d', 'temp-60d', and 'temp-90d', respectively). Five parameters of precipitation and five parameters of sunlight were examined: sum of previous 1, 3, 6, 24, and 72 h of precipitation ('precip-1h', 'precip-3h', 'precip-6h', 'precip-24h', and 'precip-72h', respectively) and sunshine duration ('sunlight-1h', 'sunlight-3h', 'sunlight-6h', 'sunlight-24h', and 'sunlight-72h', respectively). The day length on the sampling day at the three sites was calculated based on the time of sunrise and sunset obtained from the Ephemeris Computation Office NAOJ website (https://eco.mtk.nao.ac.jp/koyomi/index.html.en). The appropriate number of components was determined based on the explanative performance of the model (R2Y) and predictive performance of the model by cross validation (Q2Y) using the R package 'ropls' [28]. Correlation loading plots of the environmental parameters and the targets were created using the appropriate number of components. In addition, Pearson correlation coefficients were calculated for the normalized intensity values of the 2,505 differentially expressed targets and the meteorological data, using data from all three sites (r < −0.7, 0.7 < r, P-value < 0.001).
GO and Pfam enrichment analyses were carried out in order to determine whether the different target lists were enriched for specific GO terms or protein domains. The targets of resulting lists were categorized based on 'GO biological process complete' and compared to 'Arabidopsis thaliana (all genes in database)' using the PANTHER Overrepresentation Test (released 20220202; http://www.pantherdb.org/). Genes related to cold and UV response were extracted from among the differentially expressed targets according to the GO terms 'response to cold' and 'response to UV'. Prior to the Pfam domain enrichment analysis, the single best ORF of each target was translated into the corresponding amino acid sequence using the tool 'TransDecoder', and Pfam domains of each target were identified using the tool 'rpsblast' [29][30][31]. The enrichment P-value for each Pfam domain was calculated using Fisher's exact test in comparison to the 13,385 targets with a raw signal value exceeding 100 in each pair of samples under at least one of the 39 conditions examined (number of targets > 5, P-value < 0.01).
To further investigate the focused adaptation mechanism, genes related to hormones listed in the 'plant hormone signal transduction' pathways of A. thaliana in the Kyoto Encyclopedia of Genes and Genomes (KEGG; https://www.genome.jp/kegg/) and genes related to hormone biosynthesis listed in 'hormone metabolic pathways and genes in Arabidopsis' of the RIKEN Plant Hormone Research Network (http://hormones.psc.riken.jp/pathway_hormones.html) were extracted from among the 2,505 significant targets. Genes listed in the 'starch and sucrose metabolism' pathway in KEGG and a previous report regarding Arabidopsis [32] were extracted, assuming that those genes are related to starch and sugar metabolism. Similarly, genes related to amino acids were extracted according to the genes listed in the 'biosynthesis of amino acids' pathway in KEGG.

Quantitative RT-PCR
To evaluate the reliability of the microarray data, five genes (peroxidase superfamily protein, lipoxygenase5, euonymus lectin S3, cytochrome P450 family 718, and glycosyl hydrolase superfamily protein) that exhibited different seasonal expression patterns in the microarray analysis were analyzed by quantitative RT-PCR (qRT-PCR) using all 78 samples. Primer pairs were designed based on next-generation sequencing isotigs using Oligo version 7.6 (Molecular Biology Insights, Inc., Vondelpark Colorado Springs, CO, USA) and Primer Express version 3.0.1 (Life Technologies) (S3 Table). The first-strand cDNA was synthesized from 500 ng of total RNA using a High Capacity RNA-to-cDNA kit (Life Technologies, Carlsbad, CA, USA). qRT-PCR was performed with Power SYBR Green PCR Master Mix (Life Technologies) and a StepOnePlus Real-time PCR system (Life Technologies), as described in the manufacturer's instructions. Six microliters of cDNA diluted 1/24 with sterilized water was used in a reaction volume of 20 μl per well. Melting curve analysis was performed from 60 to 95˚C, with data captured every 0.3˚C to ensure amplification of a single product. Reaction efficiency was assessed using standard curves based on 4-fold dilution series of cDNA synthesized from 2,000 ng of total RNA (1 to 1/256 dilution). Each sample was tested independently and in triplicate using all primers. Transcript abundance was normalized to a putative actin 7 gene, which exhibited a constant microarray value according to the Pfaffl method [33], and the data obtained for each time point were compared with data obtained for shoots collected on February 18 or July 12, 2014, in Ibaraki. Very similar expression patterns were obtained from both the microarray and qRT-PCR analyses, suggesting that the data obtained in this study were reliable (S2 Fig).  Fig 1D). Growth in Ibaraki was more rapid than that in Yamagata during the growth period, from May to October. In Kumamoto, vigorous growth occurred from May to June but was slower in August and September. Height growth ceased in early October in Yamagata and Ibaraki, whereas it continued until late October in Kumamoto. ANOVA of annual growth rate indicated a significant difference in height increase among the three sites (p = 0.05). In Ibaraki and Kumamoto, the cuttings grew approximately 1.6-and 1.5-fold, respectively, during one growing season. The cuttings in Yamagata exhibited less height growth than those in Ibaraki and Kumamoto, with an approximately 1.3-fold growth relative to the initial height.

Annual transcriptome dynamics and site differences
PCA was performed to obtain an overview of the annual transcriptome dynamics at the three sites (Fig 2A). The annual time series data from all sites were plotted along a circle, indicating the seasonal life cycle of Japanese cedar. PC1 explained 69.7% of the total observed variation in gene expression. The PC1 score was high in summer (from June to September) and low in winter (from December to February) (Fig 2A and 2B). PC2 explained 5.9% of the total observed variation in gene expression. The differences in PC2 may be related to the differences in transient seasons between spring (February to April) and autumn (November) (Fig 2C). Between-site variability in annual transcriptome dynamics was also revealed by PCA. The PC1 score was higher in spring (March and April) and autumn (October) at Kumamoto than at Yamagata and Ibaraki, indicating that the transcriptome remained in the growth status for longer at Kumamoto (Fig 2A and 2B). The average PC1 scores in March, April, and October at Kumamoto were 52.3, 34.4, and 54.1 points higher, respectively, than the PC1 scores at Yamagata. The PC2 score at Yamagata was higher compared with Ibaraki and Kumamoto in winter (from December to March), indicating that gene expression differed between the sites, particularly during the dormant period (Fig 2A and 2C). The difference in average PC2 score between Yamagata and Kumamoto was 48.7 during the dormant period (December to March), whereas the difference was 9.1 during the growth period (April to November).
Hierarchical clustering demonstrated two large clusters, indicating a difference between the growth period (May to October) and dormant period (November to March) (Fig 3). The samples from April divided into the two clusters; the samples from Yamagata and Ibaraki were classified in the cluster with the samples from the dormant period, and the samples from Kumamoto were classified in the cluster with the samples from the growth period. In the cluster containing growth period samples, the samples from Yamagata and Ibaraki collected in October created a sub-cluster, and the samples from Kumamoto created a cluster with other samples from the growth period. In the cluster of dormant-period samples, all six samples from November created a sub-cluster within the dormant-period cluster.

Genes differentially expressed between sites
Analyzing annual expression patterns of the microarray data resulted in the identification of 1,473, 1,137, and 925 differentially expressed targets in the 'Yamagata vs. Kumamoto', 'Yamagata vs. Ibaraki', and 'Ibaraki vs. Kumamoto' comparisons, respectively (Fig 4). The comparisons revealed a total of 2,505 differentially expressed targets (S4 Table). A total of 644 differentially expressed targets were shared between the 'Yamagata vs. Kumamoto' and 'Yamagata vs. Ibaraki' comparisons, and 308 differentially expressed targets were shared between the 'Yamagata vs. Kumamoto' and 'Ibaraki vs. Kumamoto' comparisons. In contrast, only 153 differentially expressed targets were shared between the 'Yamagata vs. Ibaraki' and 'Ibaraki vs. Kumamoto' comparisons. A total of 75 differentially expressed targets were common among all three comparisons (Fig 4, Table 1).
The results of PLSR analyses examining the correlations between the 2,505 differentially expressed targets and meteorological parameters are shown in Fig 5. The correlations between the 2,505 targets (response variables) and latent vectors and the correlations between environmental factors (predictor variables) and latent vectors are presented. Three predictive components were used for the analyses (S5 Table), as the model had high R2Y and Q2X values (0.73 and 0.67, respectively). Parameters closed to targets were positively correlated with the nearby targets, and the meteorological parameters plotted on the opposite side (origin symmetry) were negatively correlated with the targets. As shown in Fig 5, parameters related to air temperature and day length were correlated with most of the analyzed targets. In contrast, few targets were correlated with parameters related to sunlight and precipitation.  To identify the most effective meteorological parameter for each target, correlations between the meteorological parameters and expression levels of the 2,505 differentially expressed targets were calculated (S4 Table). Most of the differentially expressed targets showed the highest correlations with air temperature and day length, with 1,842 targets (73.5%) and 222 targets (8.9%), respectively ( Table 2). The parameter air temperature at  sampling time ('temp-0h') correlated with 997 targets. The parameters of long-term air temperature more than 30 days before sampling ('temp-30d', 'temp-60d', and 'temp-90d') showed high correlations with 293, 104, and 35 targets, respectively. The parameters precipitation and sunlight did not show strong correlations with any of the targets. In addition, 441 targets did not show a strong correlation with any of the parameters examined in this study. GO enrichment analyses revealed that the 2,505 differentially expressed targets were enriched in genes associated with 'response to stress' and 'response to abiotic stimulus' (Table 3). A total of 571 targets corresponded to 485 Arabidopsis genes associated with 'response to stress', and 431 targets corresponded to 381 Arabidopsis genes associated with 'response to abiotic stimulus'. Among the genes related to 'response to stress', 24 targets corresponding to 22 Arabidopsis genes associated with the GO term 'response to cold' exhibited a strong correlation with air temperature parameters, including PRP31 and F-BOX PROTEIN 7 (FBP7) (S6 Table). The 11 targets corresponding to 10 Arabidopsis genes related to the GO term 'response to UV' exhibited a strong negative correlation with the parameter temperature, including ATCSA-1, root UVB sensitive 1 (RUS1), and UV repair deficient 7 (UVR7) (S7 Table). Genes associated with the GO terms 'cellular biosynthetic process' and 'fatty acid oxidation' were also enriched among the 2,505 targets. Pfam enrichment analyses indicated that six domains were enriched among the 2,505 targets, including 'Microtubule binding' and 'Polysaccharide biosynthesis protein' ( Table 4). All targets within the 'Microtubule binding' domain showed a positive correlation with air temperature (S8 Table).
The 824 targets positively correlated with the parameters of air temperature were enriched in genes associated with 'microtubule cytoskeleton organization', 'plant-type cell wall biogenesis', 'developmental growth', and 'cell division', which may be related to organization and growth (Table 5). In contrast, 1,018 targets negatively correlated with air temperature were enriched in genes associated with 'organic substance metabolic process', 'fatty acid beta-oxidation', and 'gene expression'. No statistically significant results were obtained for targets correlated with day length.
Twenty-three targets corresponding to 20 starch-and sugar-related genes were extracted from among the 2,505 differentially expressed targets (Table 7), and these targets showed variations in expression. Expression levels of the starch synthesis-related genes phosphoglucomutase (PGM), AGPase large subunit 2 (APL2), and soluble starch synthase 2 (SS2) were strongly positively correlated with the air temperature parameters 'temp-0h' and 'temp-30d'. Expression levels of the starch degradation-related genes like sex four 1 (LSF1) and beta-amylase 4 (BAM4) were negatively correlated with 'temp-0h', and expression levels of the gene glucan water dikinase 1 (GWD1) were negatively correlated with 'day length'.

Annual transcriptome dynamics vary between sites from autumn to spring but not in summer
Our 1-year period microarray data revealed dramatic changes in transcripts over the time course and varied expression patterns between the three sites (Fig 2). Interestingly, differences in the transcriptome dynamics were revealed from November to March. During the dormant period (from January to March), the transcriptome of Japanese cedar planted in the colder site (Yamagata) exhibited a higher PC2 score than the trees planted in the warmer sites (Ibaraki and Kumamoto) ( Fig 2C). As annual ring growth of Japanese cedar is positively correlated with temperature in February and March [6], differences in the transcriptome during this season may reflect growth. The transition to the growth period (from March to April) began earlier at the transcriptome level in Kumamoto as compared with Yamagata and Ibaraki, and the transition to the dormant period (from November to December) occurred earlier in Yamagata than in Ibaraki and Kumamoto (Fig 2). This result agreed with the results of hierarchical clustering (Fig 3). The April samples from Kumamoto belonged to the same cluster as growth period samples, whereas the April samples in Yamagata and Ibaraki belonged to the same cluster as the dormant period samples. Also, the October samples from Yamagata and Ibaraki divided from the other growth period samples and created a sub-cluster. These transcriptome differences during the transition period in spring and autumn indicated that the growth period in Kumamoto was longer than that in Yamagata and Ibaraki; thus, the transcriptome differences could have contributed to the longer growth period in Kumamoto (Fig 1D). PCA of the microarray data also demonstrated continuous changes in the transcriptome throughout the dormant period, from November to March, at all three sites. Although Japanese cedar does not grow during the dormant period, these changes may reflect preparations for growth in the following spring. Only small differences in the transcriptome were observed between the sites during the growth period. The transcriptome of all three sites exhibited higher PC1 scores from June to September, and the data plotted to a similar position (Fig 2B), even though the weekly mean air temperature on the sampling date at the three sites ranged from 19.3 to 27.6˚C. Once the growth period began, the Japanese cedar transcriptome may not have been significantly influenced by temperature. In an analysis of Japanese cedar seedlings, Ujino-Ihara [34] reported that heat shock treatment (45˚C, 3 h) induced approximately 3,000 differentially expressed unigenes. She also reported that heating pre-treatment reduced the number of differentially expressed genes induced. Thus, these data suggest that the transcriptome changes dramatically when the temperature exceeds a certain threshold and/or that Japanese cedar acclimatize to heat stress during the early stages of the active growth period under field conditions as the temperature gradually rises.

Genes differentially expressed among the three sites
The comparison of 'Yamagata vs. Kumamoto' revealed the largest number of significant targets (1,473 targets) among the three comparisons (Fig 4), perhaps because these sites exhibited the largest difference in temperature throughout the year (Fig 1B). The comparisons of 'Yamagata vs. Ibaraki' and 'Ibaraki vs. Kumamoto' revealed 1,137 and 925 significant targets, respectively. The lower number of significant genes could be attributed to the intermediate temperature in Ibaraki. The 75 common significant targets among the three comparisons may play an especially important role in environmental responses (Table 1). Among the 75 common targets, the expression of eIF1A showed the highest negative correlation (r = −0.95) with 'temp-30d'. The mechanisms of translation initiation are fundamentally similar among all organisms, and two eukaryotic translation initiation factors are universally conserved, including eIF1A [35][36][37]. eIF1A is an important determinant of tolerance to NaCl stress in yeast and plants, and overexpression of eIF1A is associated with a slow-growth phenotype [37]. The higher expression of eIF1A in Japanese cedar at lower temperatures may be related to the response to cold stress and growth cassation at lower temperatures. A positive correlation (r = 0.77) was observed between long-term air temperature ('temp-90d') and expression of the gene encoding SAM-Mtase superfamily protein, which is a key enzyme in phenylpropanoid, flavonoid, and other metabolic pathways of biotechnological importance [38]. A negative correlation (r = −0.72) between 'day length' and ABA1 expression was also identified ( Table 1). The gene encoding highly ABA-induced PP2C protein 3 (HAI3) also exhibited a strong negative correlation with 'day length' ( Table 6). Expression of ABA-related genes in accordance with changes in day length during xylem formation was also observed in the cambium region of Japanese  cedar [39]. These data suggest that the expression of ABA1 is involved in controlling endogenous ABA responses to day length and regulating the timing of xylem formation. A total of 2,505 targets were differentially expressed in the three comparisons (Fig 4 and S4  Table). GO enrichment analysis indicated that the 2,505 differentially expressed targets included a significantly higher proportion of genes associated with the terms 'response to stress' and 'response to abiotic stimulus' (Table 3). These genes may play an important role in environmental adaptation. PLSR analysis indicated that the parameters of air temperature and day length subsequently or directly influenced the expression of most of the 2,505 common targets (Fig 5). To estimate the effects of the environment on the expression level of each target, we also calculated Pearson correlation coefficients, and 2,065 targets (82.4%) showed a strong correlation (r > 0.70) with the meteorological parameters examined in this study (Table 2). Air temperature-related parameters were correlated with 73.6% of the 2,505 differentially expressed targets, indicating that temperature is the dominant factor controlling the different transcriptomes between sites. This result was consisted with a study in rice [40,41], which found that seasonal temperature changes are the dominant factor controlling transcriptome dynamics in the field. Interestingly, long-term air temperature (average air temperature over a 30-day period) was correlated with the expression of 17.2% of the 2,505 differentially expressed targets. These genes may be more affected by long-term temperature trends compared with short-term fluctuations in temperature, which could be more stable. This phenomenon indicates that Japanese cedar may be influenced by temperature trends over the previous 2 to 3 months.
The 441 targets that did not exhibit a strong correlation with meteorological parameters in this study may be influenced by other environmental parameters, such as soil conditions (e.g., temperature, moisture content, and nutrients), or the lack of correlation could perhaps be explained by examining other calculated parameters, such as cumulative temperature values as reported in an analysis of FLC gene expression in Arabidopsis [42]. Although we calculated the correlation to a single environmental parameter in this study, several environmental parameters may affect the expression level of a single target, as reported in rice [40]. Calculations of the effects of individual genes of interest could perhaps provide more details regarding the impact of environmental conditions on gene expression.

Growth-related genes up-regulated as temperature increased
GO enrichment analysis indicated that the 824 targets positively correlated with parameters of air temperature included a significantly higher proportion of genes related to 'microtubule cytoskeleton organization', 'cell wall organization or biogenesis', and 'growth' (Table 5). Pfam enrichment analysis demonstrated a significantly higher proportion of targets including 'microtubule binding' domain, and all nine targets associated with 'microtubule binding' domain exhibited a strong positive correlation with the parameters of air temperature (Table 4 and S8 Table). These results indicate that Godai1 may promote formation of the microtubule cytoskeleton and growth in summer, especially in warmer sites. Nineteen targets corresponding to 17 Arabidopsis genes were related to 'microtubule cytoskeleton organization', including the TON2 gene, which encodes a putative novel protein phosphatase 2A regulatory subunit essential for control of the cortical cytoskeleton in Arabidopsis [43]. Expression of TON2 was positively correlated with 'temp-60d'; thus, this gene may play an important role in controlling cytoskeleton-related genes. A total of 73 targets corresponding to 65 Arabidopsis genes were related to 'growth', including cellulose-related (cellulose synthase A catalytic subunit 1, cellulose synthase interactive 1, and cellulose synthase interactive 3) and expansin-related genes (expansin-A4). This phenomenon suggests that growth differences between sites could be influenced by the length of the growth period.

Acclimation to harsh winter conditions at the colder site
Among the 2,505 differentially expressed targets in this study, 24 targets related to cold response exhibited a strong correlation with the parameters of air temperature, including PRP31 and FBP7 (S6 Table). PRP31 is necessary for pre-mRNA splicing and regulation of the expression of cold-responsive genes [44], and FBP7 is required for protein synthesis during temperature stress in Arabidopsis [45]. These genes may contribute to cold acclimation in winter, especially in colder regions. Our data also indicated the possibility of an important role for chitinases in over-wintering, especially in colder regions. Among the 2,505 differentially expressed targets, a sequence of six targets exhibited high homology to the genes encoding chitinase A, homolog of carrot EP3-3 chitinase, class V chitinase, and chitinase family protein.
The expression of these genes exhibited a strong negative correlation with the parameters of air temperature (S4 Table). Plant chitinase has been implicated in defense responses against both biotic and abiotic stressors [46][47][48][49]. Diverse chitinases have been identified in spruce, and these enzymes act in concert to protect against freezing injury, store nitrogen, and induce growth cessation by promoting cell maturation [50]. Strong light in low temperature creates an imbalance between light energy absorption and energy use, leading to light stress in winter [51]. Eleven targets related to UV response exhibited a strong negative correlation with the parameters of air temperature, including ATCSA-1, RUS1, and UVR7 (S7 Table). In Arabidopsis, ATCSA-1 plays an important role in UV-B tolerance and genomic integrity [52], and RUS1 and UVR7 mutants are sensitive to UV light [53,54]. These genes may contribute to photoprotection in winter, especially in colder areas. The photoprotective role of rhodoxanthin in Japanese cedar during cold acclimation is well known [55][56][57][58]. Several photoprotective processes that prevent severe damage during exposure to strong light in low-temperature conditions may function in Japanese cedar.
Stress-related hormones may play important roles in acclimating to harsh winter conditions, especially in colder areas. Among the 2,505 differentially expressed targets, five targets related to ABA and four targets related to ethylene showed a strong negative correlation to the parameter of air temperature or day length, indicating up-regulation of the expression of these targets in winter (Table 6). Both ABA and ethylene are known to play roles in regulating bud dormancy [59]. ABA is a major regulator of bud dormancy and promotes starch accumulation in the dormant phase by regulating the expression of starch-and sugar-related genes in grapevine (Vitis vinifera L.) buds [59,60].

Starch-and sugar-related genes exhibit varied annual expression patterns and differences between sites
Starch breakdown not only enhances cold and frost tolerance, it also contributes energy for bud breaking and shoot growth [61][62][63][64][65][66]. A previous report indicated that sugar concentration peaks in winter (January) in every part (upper, middle, and lower layer shoots and roots) of the seedlings of Japanese cedar [15]. Frost hardiness is also increased in winter, especially in colder sites [16]. Pfam enrichment analysis indicated a significantly higher proportion of targets associated with 'Polysaccharide biosynthesis protein' among the 2,505 differentially expressed targets (Table 4). In addition, 23 targets encoding starch-and sugar-related genes were differentially expressed between sites (Table 7). Among the starch degradation-related genes described in a previous report [32], the BAM4 and LSF1 genes were negatively correlated with 'temp-0h', and GWD1 was negatively correlated with 'day length'. These genes may promote starch degradation at lower temperatures in northern-latitude sites and may contribute to the difference in frost hardiness between sites. In contrast, three target genes related to starch biosynthesis [32] exhibited a positive correlation with parameters related to air temperature (PGM, SS2, and APL2), and these genes may contribute to starch biosynthesis in summer. The varied expression patterns and between-site differences in starch-and sugar-related genes may indicate that Godai1 biosynthesizes sugars and starches that are optimal for the environmental conditions.

Lower temperature may induce energy accumulation in winter for growth in spring
In winter, especially in colder areas, Godai1 may prepare for growth in the upcoming spring by accumulating energy. Brassinosteroid, cytokinin, and auxin are well known as growthrelated hormones. A total of eight targets related to these hormones showed a strong negative correlation with the parameters of air temperature in this study (Table 6). We also observed higher expression of genes related to fatty acid beta-oxidation and amino acid synthesis in lower temperatures, which could be related to generation of energy for growth. These results indicate the importance of low temperature to growth in the coming spring in a clone of Japanese cedar.
Beta-oxidation is the main pathway of fatty acid degradation [67]. GO enrichment analysis indicated that the 1,018 targets that were negatively correlated with parameters of air temperature included a significantly higher proportion of genes related to 'fatty acid beta-oxidation' ( Table 5). This result appears to be consistent with those of previous reports of peanut (Arachis hypogaea L.) [68] and rice (Oryza sativa L.) [69], which exhibited up-regulation of genes related to fatty acid beta-oxidation under cold stress. Up-regulated expression of multiple beta-oxidation genes during germination was also reported in Tung tree (Vernicia fordii, Hemsl.) [70] and Upland cotton (Gossypium hirsutum L.) [71]. As fatty acid beta-oxidation generates energy for a growing embryo [71,72], it may be an important mechanism to store lipids in winter as an energy source for growth in the coming spring, especially in colder regions.
Amino acids play multiple physiological roles in plants. For example, amino acids function as nitrogen carriers in transport systems, precursors for important metabolites, nitrogen storage molecules, stress response molecules, and signaling molecules [73]. Mori [17] reported that seedlings of Japanese cedar express approximately 20 amino acids, with citrulline, glutamate, and proline being the predominant amino acids in this species. Citrulline, a hydroxyl radical scavenger, accumulates in response to drought stress and nitrogen status in watermelon [74][75][76]. Mori [77] also reported that citrulline is a major compound involved in nitrogen translocation from roots to the shoot via the xylem in Japanese cedar. Glutamate plays a very important role in plant growth and development and in the response and adaptation to abiotic stress [78][79][80]. Proline plays a highly beneficial role in plants and also provides the cells sufficient energy to sustain rapid growth [81][82][83]. With respect to the high level of citrulline and low level of arginine accumulation, Japanese cedar may differ from other gymnosperm tree species that accumulate high levels of arginine [17,[84][85][86]. Seasonal changes in amino acids have been observed in coniferous trees [17,87,88]. In Japanese cedar, citrulline and proline accumulate during the active growth period and then decrease significantly, whereas glutamate exhibits less-marked seasonal changes [17]. In this study, the expression of genes involved in glutamate biosynthesis (ALAAT2 and ASP2) [89], citrulline biosynthesis (OCT and peptidase M20/M25/M40 family protein), and proline biosynthesis (P5CR) were up-regulated as the temperature decreased ( Table 8). The gene encoding peptidase M20/M25/M40 family protein demonstrated the highest negative correlation with 7-day mean air temperature among the six glutamate biosynthesis genes, regardless of season and site (Fig 6). This result indicated the importance of low temperature for amino acid biosynthesis. Interestingly, the expression patterns of amino acid biosynthesis genes in this study and the accumulation dynamics of amino acids reported by Mori [17] were not in accordance. The amino acid biosynthesis genes were up-regulated in winter when the temperature decreased, and the amino acids accumulated in the coming spring. It is known that some mRNAs are transcribed during cold acclimation and stored for the subsequent cold de-acclimation [90]. The up-regulation of amino acid biosynthesis-related genes at low temperature (during the dormant period) may represent preparation for growth in the coming spring. The reason for the higher expression of these genes in colder sites remains to be elucidated, however. Godai1 planted in colder regions may accumulate higher amounts of amino acids to promote primary growth in spring. Warmer winters due to climate change could reduce the expression of amino acid-related genes and affect growth in spring. Further study of amino acids in Japanese cedar could increase understanding of the adaptation mechanism during the transition from the dormant to active growth periods.

Conclusions
We revealed changes in the annual transcriptome dynamics in a clone of Japanese cedar by analyzing rooted cuttings planted in three sites differing in terms of climate conditions. We identified a total of 2,505 differentially expressed targets among the three sites that may play important roles in enabling cuttings to adapt to local environmental conditions. The expression of 2,064 targets (82.4%) was affected by air temperature and day length, suggesting that these targets are directly or indirectly (subsequently) regulated by these environmental factors. Although short-term temperature was correlated with most of the genes, long-term temperature was also correlated with some genes, indicating that temperature trends may influence the expression of transcripts over several months. Targets that exhibited a strong negative correlation with air temperature included genes possibly related to cold tolerance and energy accumulation for the coming spring, such as amino acid biosynthesis, starch degradation, and fatty acid beta-oxidation genes. Transcriptome differences between sites observed from autumn to spring may reflect these differentially expressed targets. These results indicate the importance of adaptation to winter climate conditions in planting regions and the importance of low temperature to growth in the coming spring in a clone of Japanese cedar. Understanding how Japanese cedar adapt to various climate conditions is becoming more important due to changes in global climate in recent years. The results of this study may help elucidate the underlying biological mechanisms of environmental response in Japanese cedar.
Geographically, Japanese cedar natural forests range extensively from northern Honshu (30N, 130E) to Yakushima (40N, 140E) [91], and the genome-wide genetic diversity of natural populations and breeding core collections have been reported [92][93][94][95][96]. Tsumura et al. [93] studied the relationships between genotypes of natural populations and environmental variables and identified loci associated with the studied environmental variables. The results of the present study together with those of previous studies suggest the possibility of divergence of environmental responses in Japanese cedar. Although only the Godai1 clone was examined in the present study, other clones may show different environmental responses. Future studies examining additional clones selected from different regions could enhance our understanding of how Japanese cedar adapt to various climate conditions.